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It is well known that the electron transfer (ET) reaction in donor-bridge-acceptor systems is 
much influenced by the bridge conformational changes. However, the importance of the dynamical 
contribution to the ET rate is still under debate. In this study, we investigate the significance of 
the dynamical non-Condon effect on the ET rate in the coherent resonant tunneling regime. To 
, deal with the resonance, we generalize the time-dependent Fermi's golden rule expression for the 

■ ET rate in the T-matrix framework. The dynamical effect is thereby expressed in terms of the time 

' correlation function of the non-Hermitian T-matrix. We discuss the property of the quantum time 

^S) . correlation function of the non-Hermitian T-matrix and construct the semiclassical approximation. 

, ^ ■ Then we perform a computational study using classical equilibrium molecular dynamics simulations 

' and quantum chemical calculations. It is shown that the ET rate can be significantly enhanced by 

, the dynamical non-Condon effects. 

o : 

Cn . I- INTRODUCTION 

. Electron transfer (ET) is one of the most important elementary chemical processesii In many cases the simple 
i' picture of a direct transfer of an electron from the donor (D) to the acceptor (A) does not apply. Instead, bridge (B) 
molecules are involved. The standard expression for the ET rate is based on the Condon approximation, in which the 
dependence of the electronic coupling between D and A on the nuclear configuration is disregarded. However, molecules 
^ undergo structural fluctuations over a wide range of timescales at thermal equilibrium. If the B conformational changes 
influence the electronic coupling, the system will experience a breakdown of the Condon approximation. In the limit 
of slow electronic coupling fluctuations, the ET rate can be expressed as the ensemble average of the instantaneous ET 
^ rate over the static B configuration. This expression captures the structural averaging effect of the B conformation 
• ^ on the electronic coupling, while neglecting the coupling between nuclear and electron dynamics. The dynamical 
non-Condon effects may become significant as the electronic coupling fluctuations occur on faster timescales. 

There have been theoretical developments and computational studies elucidating the dynamical non-Condon effect 
fS i' on the ET rate in the deep (off-resonant) tunneling regime where the energy of the initial electronic state is well 
separated from the eigenenergies of the B subspacc^^'^'^'^'ii^'^iii^ Starting with the Fermi's golden rule expression for 
I [ the ET rate, the dynamical effect is expressed in terms of the time autocorrelation function of the electronic coupling. 
^1 In a quantum mechanical treatment, the electronic coupling is a Hermitian operator, and its time autocorrelation 
00 ' function satisfies detailed balance. The calculation results in Ref. show that for the models of independently 
coupled bridge modes with coupling of the tunneling electron to the bridge modes, the dynamical non-Condon effects 
can result in substantial enhancement of the ET rate in the activated inverted region. The detailed balance property 
, is reflected in the asymmetry of the ET rate with respect to the activation energy. On the other hand, the analogous 
CO ' classical time autocorrelation function is a real symmetric function of time and thus does not satisfy detailed balance. 
' Therefore, the classical approximation cannot reproduce the asymmetry of the ET rate due to the dynamical non- 
Condon effects, but rather produces the symmetric correction to the ET rate with respect to the activation energy, 
assuming the classical Marcus expression for the Franck-Condon factorii The enhancement in the activated inverted 
region is much smaller than that in the quantum mechanical treatment. As pointed out in Ref. 0, it is possible 
to restore detailed balance in the way that the Fourier transform of the classical time autocorrelation function is 
, . multiplied by a quantum correction factor. Such semiclassical approximation schemes have been widely discussed 
over the yearsffi>i^^>i^^>il^ia^^^^^ 

The computational studies have been carried out using classical equilibrium molecular dynamics (MD) simulations 
and quantum chemical (QC) calculations to evaluate the classical and semiclassical time autocorrelation functions of 
the electronic coupling ^i^i^ii^ The results show that the characteristic timescales for the electronic coupling fluctua- 
tions are too slow for the dynamical non-Condon effects to be significant. In the time domain, a comparison between 
the decay times of the electronic coupling correlation function, and the time-dependent Franck-Condon factor reveals 
the extent to which dynamical non-Condon effects influence the ET rate. The decay time of the time-dependent 
Franck-Condon factor is typically 3 — 4 fs.'^^ On the other hand, the decay time of the electronic coupling cor- 
relation function is more than a few tens of fsM-'SdS. Therefore, the electronic coupling fluctuations are negligible 
on the timescale of the decay of the time-dependent Franck-Condon factor. A significantly short decay time of the 
electronic coupling correlation function was observed only when the electron tunneling energy was artificially brought 
into resonance with bridge-centered eigenstatesj^ However, there has not ever been confirmation of this statement for 
a natural system in the resonant tunneling regime. 
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The purpose of this work is to investigate the timescale of the electronic couphng fluctuations and the possibihty of 
the dynamical non-Condon effects in the resonant tunneling regime through a computational study. We assume that 
transport is fully coherent, that is, the electron transfers from D to A in a single quantum mechanical process. In 
the resonant tunneling regime, another contribution comes from sequential tunneling where an electron first tunnels 
into B and then, after losing its phase memory, tunnels out of B. The dynamical non-Condon effects can substantially 
increase the contribution from the coherent tunneling, particularly if the B and A electronic states are off resonance. 
We will consider such situations. 

The choice of materials are largely motivated by the previous work for solvent-mediated ET, where B is the solvenls^. 
Instead of the C-clamp molecule, we use unlinked D-A systems. In fact, the curved saturated bridge in the C-clamp 
molecule plays an important role in keeping the distance between D and A. Instead, in this work the D and A 
molecules are held fixed during the simulation. As an example of the resonant tunneling we consider the ET from the 
LUMO of naphthalene to the LUMO of tetracyanoethylene (TCNE) in benzonitrilc (PhCN). This system lies in the 
highly activated inverted region (experimental data available for the case of acetonitrile solvent2£). Therefore, if the 
electronic coupling fluctuations occur on timescale as short as the decay time of the time-dependent Franck-Condon 
factor, the dynamical non-Condon effects may cause substantial enhancement of the ET rate. 

The Fermi's golden rule expression refers to the second-order perturbation theory result for the ET rate. In the 
resonant tunneling regime, higher-order tunneling processes can come into play. In Section [IT] we generalize the 
Fermi's golden rule expression to include all high-order tunneling processes in the T-matrix framework. The time 
autocorrelation function of the electronic coupling is thereby replaced by a time correlation function of the T-matrix. 
In practice we will limit our calculations to a small part of a large system, while the effect of the surroundings on the 
relevant system will be included by broadening the energy levels of the system. In this case the T-matrix should be 
treated as non-Hermitian. In Section IlIII we discuss the properties of the quantum time correlation functions of the 
T-matrix and then describe their semiclassical approximations, which will be used to evaluate the time correlation 
functions. In Section|lV]we give the details of the computational methods. In Section|V]we first show that disregarding 
the level broadening leads to unacceptable results even in the deep tunneling regime. Then we estimate the level 
broadening for the charge separation (CS) from the LUMO of anthracene to the LUMO of TCNE and the charge 
recombination (CR) from the LUMO of TCNE to the HOMO of naphthalene in PhCN, which are in the deep tunneling 
regime, by comparing the full T-matrix results with the second-order perturbation results. Considering the similarities 
in chemical structure, we expect a similar level broadening for the above-mentioned resonant tunneling process. Taking 
this into consideration we evaluate the time correlation functions for the resonant tunneling process in Section IVII 
Then we discuss the significance of the dynamical non-Condon effect on the ET rate. In Section fVlII we conclude the 
paper. In Appendix \X\ we sketch the derivation of a level broadening term in the T-matrix formalism. 

II. T-MATRIX FORMULATION 

In the Born-Oppenheimer approximation, the system under consideration can be described by the Hamiltonian 

Hs = (-Bd+^d)|D)(D| + (SA+i?A)|A)(A|+7?B + (VAD|A)(D| +h.c.), (1) 

with 

-ffD -5I-^D„|Dz;)(Dv|, (2) 

V 

i?A -^^^A^„|A«;)(A^i;|, (3) 

W 

a 

^AD =5m^A^„^D.a|&)|Au;)(D^;|(a|. (5) 

v,a w,b 

Here |D) and |A) are the D and A electronic states, respectively, |Dw) and \Aw) are vibrational states associated with 
the D and A diabatic energy surfaces, respectively, and |a) and \b) are B vibrational states that modulate only the 
tunneling barrier between the electronic |D) and |A) states. The D and A electronic states depend parametrically 
on the nuclear coordinates. Because of the fact that the ET occurs only at a specific configuration of the former 
vibrational modes, the Condon approximation is applicable for them, assuming the local configuration regime. On 
the other hand, since the latter vibrational modes affect strongly the electronic coupling, the Condon approximation 
cannot be made. Also, Ed and Ea are the diabatic energy surface bottoms (referred to as the electronic origins) of 
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the D and A states, respectively, iSot, and Eaw are the energies of vibrational states \Dv) and jAw) measured from the 
electronic origins, Ea is the energy of vibrational state |a), and Vuya.Awb denotes effective coupling matrix elements 
between two vibronic states |D)|Du)|a) and |A)|Aw)|5). 

The Fermi's golden rule based expression for the ET rate is given by 

kAw,b^Dv,a = —\{kw,b\VATi\Dv,a)\^5{Exi+ETiy+Ea -Ea~ Eaw - Eb), (6) 

where |D'i;,a) — \a)\Dv), \Aw,b) = \b)\Aw). The Condon approximation for the \Dv) and \Aw) states allows us to 
write the effective coupling matrix elements between two vibronic states |D) |Dw, a) and | A) | Aw, b) as a product of the 
electronic coupling element (6|VAD|a) and the nuclear overlap function (Aw|Dti), that is, 

{Aw,b\VAB\T)v,a) = (6|VAD|a)(Au.|Di;). (7) 

The electronic states can be constructed by using the Hartree-Fock orbitals at fixed nuclear coordinates. The full 
orbital space is divided into subspaces: the D subspace containing the D orbitals, the A subspace containing the A 
orbitals, and the B subspace containing the B orbitals. The partitioning Fock matrix is 

Fdd Fda Fdb\ 

Fad Faa Fab , (8) 
^Fbd Fba Fbb / 

e.g.. Fad is the submatrix of F with elements , where i belongs to the subspace D, and j belongs to the subspace 
A. 

In the deep tunneling regime, within the INDO approximation that we will use, the electronic coupling element 
between the D and A states is expressed by the effective Fock matrix element defined as 

hl = {MF\M+y2(MF\M{M ^ ^ J cbb,){<Pb,\F\M, (9) 



bib J 



Etun — F 



where F is the Fock operator, -Etun is the electron tunneling energy, |(/)d) and |0a) are the D and A orbitals of 
the tunneling electron, \4>bi) denotes the intervening B molecular orbitals, and "^b refers to a summation over the B 
orbitals bi. Strictly speaking, the effective Fock matrix element should be multiplied by a constant factor dependent on 
the D and A electronic configurations. We will consider the electronic coupling only between closed-shell Hartree-Fock 
ground states and singlet spin-adapted configurations that arise as a result of single excitations from them. While 
the constant factor is 1 for the coupling between the singly excited electronic configurations, it is y/2 for between the 
excited and ground state. 

The first and second term in Eq. ^ gives the through-space and through-bridge contribution, respectively. The 
through-space coupling becomes negligible at the large DA distance discussed below. The coupling matrix element 
is very much unaffected by the tunneling energy i^tun as long as it remains in the deep tunneling regime. In this 
case the tunneling energy can be safely varied around the average value of the D and A orbitals. However, this 
ambiguity becomes a serious problem in the resonant tunneling regime. In particular, the setting, in which the energy 
denominator in Eq. ^ becomes zero, causes the divergence of the effective Fock matrix element. In the following, 
we reformulate the ET rate expression in the T-matrix framework whereby the tunneling energy is fixed at the initial 
state energy of the tunneling electron. 

The procedure starts by writing the system Hamiltonian as = Ho + V with 

/Fdd \ 
Ho - Faa , (10) 



f Fdd 











Faa 





V 





Fbb 





Fda 


Fdb\ 


Fad 





Fab 


Fbd 


Fba 


y 



V = Fad Fab , (11) 
\Fbd Fba / 

where Hg and V are the matrix representation of the operators Hq and V, respectively, and V is treated as a 
perturbation to Hq. The eigenstates of the unperturbed Hamiltonian Hq can be used to construct the D, A, and B 
orbitals. 

Defining the T-matrix as 

f = V + VGs{E,)V, (12) 
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with the retarded Green's function 



GsiE) = , (13) 

E-Hs+iO+ 

one obtains the transition rate from the initial state \4>i) to the final state \4>f) as 

kf^,^^\{cbf\m)\'6iE,~Ef). (14) 

Here the two states and \(j)f) are eigenstates of the unperturbed Hamiltonian Hq, and Ei and Ef are the eigenen- 
ergies of states and respectively. We will consider the situation that the initial and final electronic state is 
described by the D and A orbital |(/)d) and |0a), respectively. In this case Ei and Ej are the energies of the the D 
and A orbitals that will be denoted by en and eA, respectively. The transition is induced by the perturbation V that 
couples D, A, and B. 

Strictly speaking, the nuclear kinetic energy terms should be added to the system Hamiltonian Hg because the Fock 
matrix only contains a contribution from the potential surface for the nuclear motion. Then, the delta function in 
Eq. (|14p can be rewritten in the same form as that in Eq. ([6|). Assume that the nuclear kinetic energy is the same in 
the initial and intermediate states to keep the T-matrix independent of the nuclear kinetic energies. This assumption 
seems to be reasonable under the Born-Oppcnheimcr approximation. Then, keeping in mind that the T-matrix is an 
operator in the space of the B vibrational states, the transition rate can be rewritten in the form 

kAw,b^^v,a = —\{b\Tp^T,\a)Wkw\Bv)\^5[E^ + Sd„ + Ea-E^- E^n, - Et), (15) 
n 

with Tad — ('/'a|T|(/'d)- In the deep tunneling regime, expansion of the T-matrix element Tad to second order in 
the perturbation V gives the effective Fock matrix element with substitution of the D orbital energy eD for the 
tunneling energy T'tun- The T-matrix includes all higher-order tunneling processes where the initial and final states 
are coupled by multiple scatterings described by the perturbation V. 

The observed ET rate is given by summing over final vibrational states and thermal averaging over initial vibrational 
states as 

(16) 

v,a w.b 

where Pd-u and Pa are the Boltzmann distributions over the vibrational levels |Dw) and |a), respectively. We denote 
the thermal average over the vibrational states associated with the D diabatic energy surface by (. ..)d, that is, 
(• ■ ■>D = T,v Pdv{Gv\ . . . \Dv) = TrD[e-''^° . . .]/TrD[e-''^°], where /? = {keTy^ and Tro denotes a trace over the 
eigenstates of H-a, and the thermal average over the B vibrational states by (. . .)b, that is, (. . .)b = X)a Pa{a\ ■ ■ ■\a) = 
TrB[e~''^^ . . .]/TrB[e~''^^], where TVb denotes a trace over the eigenstates of Hb- Then, writing the delta function 
in Eq. (jlSp as a Fourier transform, we get 

fcA^D = ^ j depFciEo -EA-e) CT(e/fi), (17) 

with the thermally averaged Franck-Condon weighted density of states, pFC, given by the Fourier transform of the 
thermally averaged time-dependent Franck-Condon factor 

such that 

PFc(e) = ^ y. die-*/ WW, (19) 
and with the spectral density Ct given by the Fourier transform of the time correlation function 

CT{t)^{fl^{t)fABmB, (20) 

where TadI^) = eJ-^^'^/^TA'De-'"^'^/^, such that 

CT{e/n) = / dte'^'/^Crit). (21) 



5 



Applying the convolution theorem for Fourier transforms, Eq. (|17p can be written as 

1 1"°° 

kA^^^-J dte'(^--^-)*/VFc(i)C^T(i). (22) 

Note that the T-matrix is non-Hermitian due to the infinitesimal term iO~^ . In the deep tunneling regime it is treated 
as a Hermitian matrix by disregarding the iO'^ term. However, we will see that in the resonant tunneling regime it 
is necessary to modify the eigenvalues of in the Green's function (|13p to include a finite imaginary term, which 
represents the effect of the surroundings on the relevant system. Then it is important to deal with the non- Hermitian 
T-matrix. In this case the operator Tad is also non-Hermitian. 

III. TIME CORRELATION FUNCTIONS OF T-MATRIX 

The non-Hermitian operator Tad can be decomposed into the sum of two Hermitian operators as 

Tad = i + tB, (23) 

where A and B are Hermitian operators given by 

A=^{fAB+fl^), (24) 

S=^(TAD-TiD)- (25) 

Then the quantum time correlation function Cxit) can be separated into the autocorrelation and the cross-correlation 
part, that is. 



with 



where 



CTit)^CA{t)+iCcit), (26) 



CA{t) = CAA{t) + CBB{t), (27) 

Cc{t) ^ C'ABit) ~ CsAit), (28) 

CAA{t) = {A{t)A{Q))B, (29) 

CBBit) = {B{t)B{0))B, (30) 

CAB{t) = {A{t)B(Q))B, (31) 

CBA{t) = {B{t)A{Q))B, (32) 

A{t)^^{fAD{t)+fUt)), (33) 

B{t) = ^^{fAu{t)-fl^{t)). (34) 



These two parts have the following properties: 



CA(-t) = C*A (t) = CA{t- ipn) , (35) 
Cci-t) = -C*c{t) = -C'cit - tl3H). (36) 



Thus, the Fourier transform 



satisfies the detailed balance relation 



oo 

iiot 



CAiuj) = / dte"^*CAit) (37) 



Ca{uj) = e^f'^CAi-u;), (38) 
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while the Fourier transform 

/oo 
dte'^'Ccit) (39) 
-OO 

docs not satisfy detailed balance but 

Cc{^)^~e'""^Cchij). (40) 

The cross-correlation part represents the nonequilibrium properties of the system. In the deep tunneling regime, 
Tad is often treated as a Hermitian operator, assuming that the D and A electronic wave functions are real, so that 
the cross-correlation part vanishes. From Eqs. ((551) — ([55]) it follows that C'a(w) and C'c(w) are real and imaginary, 
respectively. While C'aI'^) > 0, iCc{(^) can also be negative. Since the trace is invariant to cyclic permutation of the 
operators inside it, and the equilibrium density operator e~^^^ and the time evolution operator e~«^fBt/fi commute, 
the quantum time correlation functions are stationary, that is, {A{t + t') A{t'))B — {A-{t) A{Q))b, etc. 
It will be convenient to decompose CA{t) and Cc{t) into the real parts 

ci+\t)^^{CA{t) + CX{t)), (41) 

C^c^\t)^liCcit) + Cm), (42) 

and the imaginary parts 

ci-\t)^^{CA{t)~CX{t)), (43) 

C^c\t) = ^iCc{t)~Cm)- (44) 
From Eqs. ((35| — (|39| the Fourier transforms can be written as 



ci^\Lo)^-iCAiLo)±CA{-co)), (45) 



4±)(^) = i(C'c(a;)T^c(-^)). (46) 

While Cj^\lu) and ■'(tj) are symmetric in w, c'}^\lj) and c[^\lj) are antisymmetric in w. Substitution of 
Eqs. ^ and (gD]) into Eqs. (gS]) and (glD leads to 

^A(-) = Y^|r^cf^(-), (47) 
^c(-) = Y^|:^4^^(-). (48) 

The quantum time correlation function Ct (t) can be evaluated using classical equilibrium MD simulations coupled 
with QC calculations. This approach involves the replacement of the quantum time correlation function by its classical 
counterpart. We denote the classical time correlation function by 



with the classical variables 



Cf^(t) = (TlD(t) Tad(0))c1, (49) 

TAuit) = (0A(QB(t))inQB(t))|(/'D(QB(t))>, (50) 

TiDit) = {MQBitmHciBitmAiQBm- (51) 

Here, the T-matrix T(QB(t)), and the D and A orbitals |(/'D(QB(i))) and |0A(QB(i))) are parametric functions of the 
B nuclear coordinates Qb(0 at time t that move along the classical trajectories obtained by the MD simulations, and 
the classical thermal average (. . .)ci should be understood as the time average over the MD trajectory, as implied by 
the ergodic theorem of statistical mechanics. 

Using the classical counterparts of CA(i) and Cc{t), the classical correlation function C^{t) can be written as 

C^At) = Ctit)+^C^^{t), (52) 
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with 



Cf(t) = C''L(t) + C'ilr.(t), 


(53) 




('54) 


C^'^(t) = (A(t)A(0))el, 


(55) 


t^BBW = (SWS(0))cl, 


(56) 




(57) 




(58) 



where 

A{t) = ^{TAu{t)+n^{t)), (59) 

i?(t) = :^(TAD(t)-TlD(i)). (60) 

From the time reversal symmetry of the classical equations of motion, it follows that 

Ct{t)=C^^i-t), (61) 
C£'(i)=CcM)- (62) 

Furthermore, the stationary property implies that 

= -ClKt). (63) 

Alternatively, taking the classical limit h ^ 0, Eq. ([36]) yields Eq. (f63|) . From Eqs. ([62|) and ((63|) it follows that 

C£'(t) = 0. (64) 

Namely, the classical approximation disregards the contribution of the cross-correlation part Cq (t) in the time corre- 
lation function CT{t). 

Using Eqs. (1521), §M, dSS]), ([55]), and §UI yields 

Cf?(i) = Ctit) = i(rAD(t)rlD(0) +TlD(i)TAD(0))cl 

= (ReTAD(ORcTAD(0))ci + (InirAD(t) ImTAD(0))ci 

= QeW+QtW, (65) 

where Crl^{t) and C£!jj(t) define the autocorrelation functions of the real and imaginary parts of Tad, respectively. 
Since C'^{t) is a real symmetric function of t, its Fourier transform 

/oo 
dte'"*Cl'(t) (66) 
-oo 

is a real symmetric function of uj, that is, C'^{lu) — C'^{—u!). This means that the classical time correlation function 
does not satisfy the detailed balance condition ([55)1 . The simplest way to solve this problem is to identify C'^{t) with 

C'^\t). Using the relation (gT]) we obtain 

C'a(c^) = Q(c^)Cl'(^), (67) 
where Q{uj) is the quantum correction factor given by 

which guarantees detailed balance. There are other approximations involving quantum correction factors such as 

Qi-) = T^k' (69) 

Qiuj) = e^"^/^. (70) 
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More complex examples can be found in Refs. |21||22||23||24| 

Employing such a semiclassical approximation, the ET rate P?]) can be expressed as 

kA^u^^J depFc{EB-EA-e)Q{e/h)C'^{e/h). (71) 

Here, from Eqs. ((65| and ((66|) . C'^{uj) = C^^{uj) + C'f^(w). The general expressions for PFci^) have been widely 
discussed.— 1^ Adopting the simplest classical expression 



where A is the reorganization energy, Eq. (|7ip yields 



PFc(e) = ^/^-T^exp -^-^^ ] , (72) 



-D 



which reduces to the form given in Ref. 10 if we replace Tad by the real terms up to second order in perturbation 
theory. According to Eq. (j72p . pfc{^) has a maximum for e = A, and the energy scale for the decay is given 
by Aepc — \/4AfcsT. If C'^{e/h) goes asymptotically to zero for |e| larger than the width Aca ^ Aepc, then 
PFC {Eb — Ea — e) in Eq. ([73]) can be replaced by the value at e = to yield the static limit expression 

= ^Pfc{Eb-Ea) J deQ{elh)Cf{elh). (74) 

The corresponding quantum expression is given by 

fcA^D = -j-f>Fc{EB ~ Ea) Ct(0). (75) 

Taking this limit corresponds to replacing 5{Efi+ Ef)v + Ea — Ea — Eaw — Eb) in Eq. ([TS]) by 6{Efi +Ef)v — Ea — Eaw)- 
In this case the initial and final vibrational energies of B are the same. Hence, in Eq. (|75p Ct(0) should be replaced 
by J2a ^1 I (^I^adIa) P^Bb,^^ • When Aca is comparable to Aepc, the dynamical corrections become significant. As 
shown in Eqs. ([55)1 — (TfO)) . the quantum correction factors Q{uj) are positive increasing functions. Moreover, C'^{uj) > 0. 
Thus, due to the convolution in Eq. ([73]), the dynamical non-Condon effects can result in substantial enhancement 
of the ET rate in the activated inverted region. If the dynamical non-Condon effects are significant, its resultant 
ET rate is highly asymmetric with respect to the activation energy. In Eq. ((73)) . Q{e/h) C'^{e/h) is the semiclassical 
approximation of the spectral density CTi^/fi), which represents the electronic transition rate between D and A, 
averaged over the thermal distribution of B vibrational states, and e denotes the energy difference between the D and 
A diabatic vibronic states |D)|D?;) and |A)|Az«), that is, e = Ef, + Ef,v — Ea — Eaw, which is equal to the energy 
transferred to/from B during the electronic transition. This can be understood by rewriting pfc{Ef> — Ea — e) and 
Crie/h) in Eq. ^ as 

PFc{Eu -EA-e)^J2 Pdv\{Aw\Dv)\^6{Eu + Ej,, - Ea - Ea^ - e), (76) 

V,W 

CTie/H) = 2T:hJ2Pa\{b\TAB\a)\^S{Ea~Eb + e). (77) 

According to the Franck-Condon factor, downward (i.e. e > 0) electronic transitions are more probable than upward 
(i.e. e < 0) electronic transitions in the inverted region, while the opposite is true in the normal region. The 
asymmetry of the ET rate reflects the detailed balance property that downward transitions are more probable than 
upward transitions. This is a feature of the autocorrelation part CA{t) of Crit). If the cross-correlation part Cc(i) 
contributes to the ET rate, the dynamical non-Condon effects can result in a more complicated dependence of the 
ET rate on the energy gap Efi — Ea due to the minus sign in Eq. ([^0]) . which docs not exist in the detailed balance 
relation. In any case, the dynamical non-Condon effects are accompanied by the vibrational excitation or deexcitation 
of B and in this sense inelastic. 

Within the semiclassical approximation Eq. (I7ip can be formally written as 

1 /"OO /"OO 

fcA^D^Ti/ die'(^--^^)*/VFc(0 / dt'Q{t')Ct{t-t'), (78) 
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with Q{t) the inverse Fourier transform of the quantum correction factor, defined by 

Q{t) = — / dwe-*"*g(w). (79) 

Performing the inverse Fourier transform of Eq. (j72p . we obtain 

ppcit) = e-^'=-^*'/'^'e-'^*/^ (80) 

From this the characteristic decay timescale is given by 

TFC = h/^XkBT. (81) 

If Tpc is much smaller than the shortest timescale ta on which C'^{t) changes, we can replace C'^{t — t') in Eq. ((78|) 
by CX'(i') to get the expression (I74p . When ta is comparable to rpc, the dynamical non-Condon effects will become 
significant. In the following, we attempt to evaluate C'^{t) in the resonant tunneling regime, using MD simulations 
and QC semiempirical calculations. 

IV. COMPUTATIONAL METHODS 

We consider the ET processes between anthracene and TCNE and between naphthalene and TCNE in PhCN 
solvent, where the solvent molecules act as B. 

The MD simulations were carried out by using the program NAMD.-- We also use the program WORDOM to 
analyze the MD trajectoriesi^ We employed the AMBER force field, '^^''^^ converted for use in CHARMM, and used 
the Grabuleda et al. force field^ for acetonitrile in place of that for the cyano group of PhCN. To determine atomic 
charges of the D, A, and solvent molecules, B3LYP/6-31G(d) optimizations of individual molecules were performed, 
and then the ESP fits were implemented using Gaussian 03^. The temperature was maintained at 300 K by means 
of Langevin dynamics. Also, the particle-mesh Ewald method for full-system periodic electrostatics was taken into 
account. The D and A molecules were placed in the middle of a solvent cube and face to face with a 7 A distance in 
such a manner that the C(9)— C(10) line of anthracene or the central C— C single bond of naphthalene and the central 
C=C double bond of TCNE were parallel to each other, after energy minimization. The number of solvent molecules 
was 550. To build the starting geometry for the production dynamics, the system was equilibrated at 300 K and 
1 atm with periodic boundary conditions. The system was first energy minimized, and then a 100 ps equilibration 
dynamics led to a starting simulation box of c.a. 46.2 x 46.2 x 46.2 A. The simulation was performed in the NVT 
ensemble where the number of particles, N, the volume, V, and the temperature, T, are kept constant. Hydrogen 
atoms in the solvent molecules were constrained using the SHAKE algorithm, and also the D and A molecules were 
held fixed during the simulation. The remaining degrees of freedom were left flexible. The integration time step was 
1 fs. The length of each trajectory was in the 50 ps range, and the interval between QC calculations was in 1 fs range. 

The PM3 methoct^^ in the M0PAC7 program was used to perform the QC calculation. For each snapshot of the MD 
trajectory, the atomic coordinates of the D, A, and solvent fragments are extracted and used for the QC calculation. 
The solvent molecules involved in the coupling was selected on a geometrical basis as the ones entirely within or 
partially within a 3.5 A radius from the center of mass of the D and A molecules. 

V. LEVEL BROADENING 

We limit our calculations to the relevant subspace of an overall system. Then we need to characterize the effect of 
the rest of the system. This effect can be included by introducing self-energy terms in the denominator of the Green's 
function (fT3|) for the isolated relevant system. Considering the weak perturbation V, we expect that the second-order 
perturbation theory provides sufhcient accuracy in the deep tunneling regime where the energy of the initial electronic 
state is well separated from the eigenenergies of the B subspace. Then the self-energy terms can be safely disregarded. 

In Figs. [Hand [3 we show the second-order perturbation results for (a) the CS from the LUMO of anthracene to the 
LUMO of TCNE, and (b) the CR from the LUMO of TCNE to the HOMO of naphthalene, which are in the deep 
tunneling regime. 

Figures[TJa) and (b) show the patterns of Tad(0 second-order perturbation theory, denoted T^Q*(t), for portions 
of the MD trajectories with anthracene and naphthalene, respectively, and the corresponding time correlation functions 
Ca(Oi normalized to their initial values, are plotted in Figs.[2Ia) and (b). We denote the normalized time correlation 
function by Cf{t) = CX'(i)/Ci'(0). 
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The effective correlation time ta niay be taken as the time for which the function C'^{t) — C'^'(oo) decays to 1/e of 
its initial value. From Figs.[H^a) and (b) we can see that ta is more than a few tens of femtoseconds in agreement with 
the previous observations . Using the value A = 0.902 eV, from Table 2 of Ref. [H, in Eq. ([8T|) yields tfc ~ 4.3 fs. 
Thus, TA is too large to introduce significant dynamical corrections to the ET rate for these processes. 

Figures[3l^a) and (b) show the results of full T-matrix calculations on the same systems. Comparison of Figs. [l]and 
[3] indicates that the full T-matrix calculations gives a significant error, which mainly results from the small difference 
between the initial electronic state energy Ei and the eigenenergy of the Hamiltonian Hg in the Green's function 
(I13p . This error can be corrected by introducing the self-energy terms, which describe the effect of the surroundings 
on the relevant system. Thereby the free Green's function Gg (E) in the T-matrix is replaced by the effective reduced 
Green's function (see Appendix . 

The self-energy is a non-Hermitian matrix. The anti-Hermitian component of the self-energy is responsible for the 
broadening of the energy levels, while the Hermitian component can conceptually be viewed as a correction to the 
system Hamiltonian Hs ■ For simplicity we ignore the mixing of the states of the relevant system by the coupling to the 
surroundings. In this case the self-energy becomes diagonal in the basis of eigenstates of the system Hamiltonian Hs- 
Furthermore, the energy dependence is disregarded, and all the diagonal elements are set equal. Consequently, the 
real part of the self-energy vanishes, and the reservoir effect is characterized by a single level broadening parameter, 
which we denote F. Due to the level broadening term. Tad is complex, and the time correlation function C'^{t) is 
decomposed into two autocorrelation functions of the real and imaginary parts of Tad as in Eq. (|65p . Figures |3] and 
[5] show the consequences of including the level broadening. 

Figure 2] summarizes a representative set of simulations of RcTad (t) (green line) and ImTAD (t) (red line) for the 
CS from the LUMO of anthracene to the LUMO of TCNE. For (a), (b), (c), and (d), the level broadening is chosen 
as (a) F = 1 eV, (b) F = 0.1 eV, (c) F = 0.01 eV, (d) F = 0.001 eV. For comparison Tl^\t) of Fig.[T]is superimposed 
on each graph (black line). 

Reasonably good agreement of the full T-matrix result with T^D*(i) is obtained by choosing F to be ~ 0.01 — 0.1 
eV (Figs.m^b) and (c)), where the contribution of ImTAolO is negligibly small compared to that of ReTAD(0- When 
F is set larger or smaller (Figs.|4l[a) and (d)), the contribution of ImTAD(0 becomes significant, and the full T-matrix 
result apparently deviates from T^'^^{t). 

We get similar behavior patterns for the CR from the LUMO of TCNE to the HOMO of naphthalene as shown in 
Fig. [51 though the deviation at small level broadening becomes smaller due to the higher energy separation between 
the D and B LUMO levels. 

To analyze the origin of the deviation, we compare the full T-matrix and the second-order perturbation results for 
the CS from the LUMO of anthracene to the LUMO of TCNE, including the level broadening in both calculations, in 
Fig. [HI The full T-matrix includes all possible intermediate multiple-scattering processes and thus the D and A orbitals 
as intermediate states. On the other hand, the second-order perturbation theory exclude them from the intermediate 
processes. Consequently, the deviation of the full T-matrix result from the second-order perturbation result arises due 
to the participation of the D and A orbitals, particularly the initial electronic state, in the intermediate processes. In 
Fig. [HI the gray line is the second-order perturbation result, and the black line is the full T-matrix result. Figures [Hljc) 
and (d) clearly show that the deviation becomes apparent when we use a smaller broadening parameter. In Figs. [IK a) 
and (b) both results are almost identical. Similar results can be obtained for the CR from the LUMO of TCNE to the 
HOMO of naphthalene (data not shown). In conclusion, the deviation between the full T-matrix and the second-order 
perturbation results at small level broadening can be attributed to an underestimate of the level broadening for the 
initial electronic state, while the deviation at large level broadening can be attributed to an overestimate of the level 
broadening for the B states. In the latter case significant differences do not appear between the full T-matrix and 
the second-order perturbation calculations including the level broadening because they both involve the B states for 
which the large level broadening causes the error. 

VI. RESONANT TUNNELING REGIME 

We now turn to the resonant tunneling regime. In the previous section we have considered the CS from the LUMO 
of anthracene to the LUMO of TCNE and the CR from the LUMO of TCNE to the HOMO of naphthalene in PhCN 
solvent. Here the LUMO energies of anthracene and TCNE and the HOMO energy of naphthalene lie deep inside 
the solvent HOMO-LUMO gap. To realize the resonant tunneling regime, we consider the CS from the LUMO of 
naphthalene to the LUMO of TCNE in the same solvent, whereby the LUMO energy of D is shifted into resonance 
with the solvent LUMO. In the previous section we have seen that for F ~ 0.01 — 0.1 eV we get the best results. This 
estimation of the level broadening may also be used for the D, A, and B states here, considering the similarities in 
chemical structure, namely, cyano groups in TCNE and PhCN, and aromatic rings in anthracene, naphthalene, and 
PhCN. These similarities suggest similar interactions with the reservoir, and thus we expect a similar level broadening. 
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Figure [7] shows the patterns of Re Tad (i) (black line) and ImTAriit) (gray line) with the level broadening parameter 
chosen as above in (a), (b), (c), and (d), respectively. From this we can see that the T-matrix includes a substantial 
contribution from the imaginary part. In Fig. [8] the same data (black line) are shown in comparison with the second- 
order perturbation results with the level broadening (gray line). This indicates that for the small level broadening 
the second-order perturbation theory fails badly in the resonant tunneling regime. In the following, we consider only 
the full T-matrix. The black solid lines in Fig. [9] show the results for the autocorrelation functions C^^{t) and Cf^{t), 
normalized to their initial values. We denote the normalized autocorrelation functions by C^^{t) = CyIcW / C^^{0) 
and Cfl^{t) — Cj^(i)/Cf^(0). The level broadening parameter in (a)— (d) is chosen as above. For (b), (c), and (d), 
^Rc(*) exhibits very fast initial decay and subsequent periodic damping with a period of approximately 20 fs. In all 
cases, the oscillation frequencies are almost equal, while the amplitudes of C^^{t) are larger than those of Cf^{t). 
This feature can be understood by analyzing the dominant contributions to C^{t) and C^^{t). The gray solid lines 
in Fig. [HI show that for C^g{t) and C'f^(t) in (b) and (c) and C'f,'„(i) in (d) the dominant contribution comes from the 
two adjacent eigenstates of the system Hamiltonian Hs, which are closest to the LUMO of D (naphthalene), namely 
the initial electronic state and the LUMO of B (PhCN solvent). The deviations in (a) can be explained by the 
large contributions from the other eigenstates due to the large level broadening. On the other hand, the deviation 
from C^^{t) in (d) is more subtle. We will return to this issue later on. 

We denote the two adjacent eigenstates by and with the higher eigenenergy and the lower eigenenergy 
E-, respectively, their contribution to TAi){t) by T^|)(i), and the contributions to C^^^{t) and Cf^{t) by C^^^{t) and 
Cf^'^{t), respectively. The two states and jV'-) alternately fluctuate back and forth between the LUMOs of D and 
B, which are nearly degenerate. To demonstrate this more clearly, we consider a two-state system whose Hamiltonian 
is written as a sum Hs = Hq^ + V. Here we assume that the eigenfunctions of Hq^ are given by the initial electronic 
state \(pi) and the solvent LUMO, denoted by \(I>b), with the corresponding eigenvalues Ei and eb, respectively. These 
states are coupled by the perturbation V. For simplicity we assume that the two-state system is unaffected by the 
coupling between the B and A orbitals \(f>B) and |0a), although the perturbation V includes its contribution. In this 
case, in the basis of the functions \(f>i) and |0b), ^^s is represented by the matrix 

*s = (,5, "1 , (B2) 



Ve''' eB 

where we have denoted Ue*'' = (0B|V^|(^i) with V taken real and positive and the phase factor e"'. The eigenstates 

and eigenvalues of the Hamiltonian Hs are denoted \tp+) and \tp-), and £'+ and E_, respectively. The eigenvalues 
are given by 

E, . (83) 

Thus, \ip+) remains above the initial electronic state and the solvent LUMO while \ip-) remains below them. The 
eigenstates are expressed by 

^ cose e-'''/^\(l>^)+ sine e'''/'^\(j)B), (84) 

I^A-) = -sin6le-"'/V.:> +cos6'e"'/VB), (85) 



where 



sin<? = 7^^^' (86) 



(l+X2)l/2^ 
1 

(1+^2)1/2^ 



cose = 7^^^' (87) 



^_E+-E, V V 



V E+-eB E, - E^ 

Here V and the eigenstates and eigenvalues of Hq'^ and Hs depend on time via the solvent nuclear coordinates. 
Figure [TO] shows that E+ (dark gray solid line) and E- (gray solid line) are very good approximations to the exact 
eigenenergies E+ (black dashed line) and E- (black dotted hue), respectively, and Fig. [TT] shows that sin^ (dark 
gray solid line) and cos^ (gray solid line) reproduce the population exchange of \(f)i) and |(/)b) between the exact 
eigenstates and \ip-), described by (black dashed hue) and (black dotted line) or |((/)b|V'+)P 
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(black dashed line) and |(0b|'0-)I (black dotted line). These results allow us to approximate by 

1 



Ti2s 
AD 



E, -E++iT/2 
1 



E,-E. 



zr/2 



(V^+I^I0d) 



Since by definition (0i|y|(/)D) = 0, we have 



(89) 



(90) 



with 



(£;, - S+) sin^ 6* (E,-E^) cos 



{E,-E+f + {T/2Y {E,-E_Y + {Tl2f 



{Mv\(j)B){MvW 



(91) 



and with 



ImT, 



AD 



(r/2) sin^ d 



(r/2)cos2 6l 



.{E, - E+Y + (r/2)2 [E, -E^Y + (r/2)2 

/+ + /-. 



{c^f,\v\<PB){Mv\M 



(92) 



Here, Rj^ and i?_ denote the contributions of the states 1-0+) and respectively, to the real part ReT^p, and 

and /_ denote the contributions of the states and respectively, to the imaginary part ImT^^. Since 

E+ > Ei > E-, R+ and i?_ have opposite signs while /+ and /_ have the same sign. The black dashed lines in 
Fig. [9] show this approximation results for C^^^{t) and Cf^^{t). These plots show very good agreement with the exact 
results (gray solid lines) except for (d). The deviations are exposed due to the small level broadening. According 
to Eqs. ([SS]) . ((9T|) . and ([M]) . 1-0+) and which are nearly identical to 10+) and \tp-), respectively, contribute to 

ReT^p destructively, while they contribute to ImT^I) constructively. Thus, ReTAD is more influenced by the other 
states than ImTAD- This is the reason why the deviation of C^^^{t) from C^^{t) is much larger than that of Cf^^{t) 
from Ctit). 

For the small level broadenings F ~ 0.001 — 0.1 cV, the timescale for the initial decay of C^^{t) reaches several 
femtoseconds and thus may be fast enough to introduce significant dynamical corrections to the ET rate expressions. 
On the other hand, if the level broadening is as large as 1 eV, the dynamical non-Condon effects will be insignificant. 
The rapid initial decay of C^^{t) is followed by the periodic behavior with much the same period. The fluctuations 
of RcTad and ImTAo are caused by the solvent motions. From the timescale viewpoint, we infer that the periodic 
behaviors of the autocorrelation functions are caused by bond stretching vibrations of PhCN. The higher frequency 
components in the initial decay should be created by the coupling between various vibrational modes. 

The ET rate can be calculated using Eq. (|75|) . For C^{e/ti)/C'^{0) we use the Fourier transforms of C^^{t) and 
Cf^{t) in Fig. M Figure [TO] shows the energy gap dependence of ^a^d for the level broadening chosen as above 
in (a), (b), (c), and (d), respectively. The gray solid, black solid, and black dashed lines are the results for the 
quantum correction factors of Eqs. ([68]) . ([69]) . and ([70]) . respectively. As expected, for the small level broadenings 
F ~ 0.001 — 0.1 eV, the dynamical non-Condon effects cause the substantial enhancement of the ET rate in the 
activated inverted region. This enhancement becomes larger with decreasing level broadening. It also depends on the 
quantum correction factor used. The quantum correction factor of Eq. I|69p generates larger enhancement than that 
of Eq. ((S5|). On the other hand, that of Eq. (1701) causes extraordinarily large enhancement of the ET rate regardless 
of the level broadening because this approximation badly overestimates the value of Ca(i^) for high frequencies. 



VII. CONCLUSION 



We have investigated the significance of the dynamical non-Condon effect on the ET rate in the coherent resonant 
tunneling regime. The time-dependent Fermi's golden rule expression for the ET rate has been generalized to deal with 
the resonance in the T-matrix framework. The dynamical effect is thereby expressed in terms of the time correlation 
function of the non-Hermitian T-matrix. This time correlation function can be separated into the autocorrelation 
and the cross-correlation part. The former satisfies detailed balance, whereas the latter does not satisfy it. When the 
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shortest timescale for the change of the time correlation function is comparable to the decay time of the time-dependent 
Franck-Condon factor, the dynamical non-Condon effects become significant. 

In the semiclassical approximation the Fourier transform of the quantum time correlation function is replaced by the 
classical time correlation function multiplied by the quantum correction factor. In this case the cross-correlation part 
is disregarded. The classical time correlation function has been evaluated using the combined MD/QC simulations. 
As an example of the resonant tunneling we have considered the ET from the LUMO of naphthalene to the LUMO 
of TCNE in PhCN. For simplicity the D and A molecules were held fixed during the simulations. Experimentally, 
similar situations can be created by linked D-A systems such as C-clamp molecules. 

We have introduced the level broadening parameter to include the effect of the surroundings on the relevant system. 
Disregarding the level broadening leads to unacceptable results even in the deep tunneling regime. The parameter 
values have been estimated for the CS from the LUMO of anthracene to the LUMO of TCNE and the CR from the 
LUMO of TCNE to the HOMO of naphthalene in PhCN, which are in the deep tunneling regime. Considering the 
similarities in chemical structure, we expect a similar level broadening for the above-mentioned resonant tunneling 
system. By comparing the full T-matrix results with the second-order perturbation results, the best parameter values 
are found to be F ~ 0.01 — 0.1 eV. For such level broadenings, the timescale for the initial decay of the time correlation 
function reaches several femtoseconds, and the ET rate can be significantly enhanced by the dynamical non-Condon 
effects. 



APPENDIX A 

We consider a Hamiltonian written as a sum 

H^Hs + Hr + Vsr, (A1) 

where Hs is the Hamiltonian of the relevant system, given by the Fock matrix ([8]), 7?r is the Hamiltonian of the 
surroundings, and Vsr is the interaction between them. Let {|s)} and {\r)} be the orthonormal sets of eigenstates of 
Hs and Hr, respectively. Then Hs, Hr, and Vsr may be written in the representation as 

i7s = E£;.k)(5|, (A2) 

S 

Hn = J2Er\r){rl (A3) 

r 

Vsr = Y.^Vsr\s){r\ + (A4) 
where we denote Vsr — (^IVsrI^), etc. The retarded Green's function of the total Hamiltonian is given by 

E - H + tO+ 

We focus on the relevant subspace of the total system and attempt to characterize the effect of the surroundings on 
the relevant system. The free retarded Green's function of the relevant system is given by Eq. (fT3|) . Then the Dyson 
equation is expressed as 

G{E) = Gs{E) + Gs{E) Vsr G{E). (A6) 
Using the resolution of the identity operator 

s r 

leads to 

Gs's{E) = {Gs)ss{E)5s'^s + {.Gs)s's'{E)Y,ys'rGrs{E), (A8) 

r 

Grs{E) = {Gs)rr{E) ^ Vr G , ,{E) , (A9) 



14 



where Gs's{E) = {s'\GiE)\s), {Gs)ss{E) = {s\Gs{E)\s), etc. Inserting Eq. ^M\l into Eq. jM]) it follows that 

J2 [(E - Es' + i0+)6s',s" - Bs's"iE)] Gs"s{E) = 5s', s, (AlO) 

s" 

where Bsis"{E) are the self-energy matrix elements given by 

Bs's"{E)^y ^ . (All) 

' ' ^ ' ^ E - Er + iO+ ^ ' 

r 

This defines the effective reduced Green's function for the relevant subspace. Note that the self-energy matrix is 
non-Hermitian due to the infinitesimal term zO"*". For simplicity we assume that the mixing of the states |s) by the 
coupling to the reservoir can be disregarded, namely we employ the diagonal approximation to the self-energy matrix. 
Thereby the self-energy matrix elements read 

Bss'{E)^ Bs{E)5s,s', (A12) 

with 

Bs{E)=y (A13) 

r 

and thus the reduced Green's function becomes a diagonal matrix, that is, 

= iE-Es + !o^)-BsiE) '^'^- (^^4) 

The imaginary part of the self-energy is equal to the broadening of the density of states in the relevant system, while 
the real part corresponds to the shifts in the energy levels Eg. 

Furthermore, we assume that the set of states |r) constitutes a continuum of states. In this case the summation 
over r can be replaced by the integral 



/ dErPniEr), (A15) 
where pn{E) denotes the reservoir density of states. Equation (jA13p now takes the form 



where (|l^srp)_E is the average of the squared coupling over all continuum levels r that have energy E, defined by 

W^A^ ^Y.r\ysr?5{E - Er)/Y.rS{E - Er) = {pK{E)r^Y.r\ysr?5{E-Er), and whcrC 



Ts{E) = 2n{\VsrV)EPK{E). (A17) 
In the simplest case where Ts{E) does not depend on E and s, we get 

In this case the effective reduced Green's function can be expressed in the form 

G(E) = . (A19) 

E-Hs + iT/2 

Note that the infinitesimal term iO+ in Eq. (|A14|) can be disregarded relative to ir/2. 
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FIG. 1: Tl^\t) for (a) the CS from the LUMO of anthracene to the LUMO of TONE in PhCN, and (b) the CR from the 
LUMO of TONE to the HOMO of naphthalene in PhCN. 
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FIG. 3: TadH) without level broadening for (a) the CS from the LUMO of anthracene to the LUMO of TONE in PhCN, and 
(b) the CR from the LUMO of TONE to the HOMO of naphthalene in PhCN. 
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FIG. 4: KeTAi>(t) (green line), ImrAD(i) (red line) with level broadening, and T^'^^{t) (black line) without level broadening 
for the CS from the LUMO of anthracene to the LUMO of TONE in PhCN. The level broadening is chosen as (a) T = 1 eV, 
(b) r = 0.1 eV, (c) r = 0.01 eV, (d) T = 0.001 eV. 
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FIG. 5: RerAD(t) (green line), ImrAD(i) (red line) with level broadening, and Tj^^{t) (black line) without level broadening 
for the CR from the LUMO of TONE to the HOMO of naphthalene in PhCN. The level broadening is chosen as (a) T = 1 eV, 
(b) r = 0.1 eV, (c) r = 0.01 eV, (d) V = 0.001 eV. 
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FIG. 6: ReTAD(t) and ImTAD(i) with level broadening for the OS from the LUMO of anthracene to the LUMO of TONE in 
PhCN. The gray line is the second-order perturbation result, and the black line is the full T-matrix result. The level broadening 
is chosen as (a) F = 1 eV, (b) T = 0.1 eV, (c) T = 0.01 eV, (d) T = 0.001 eV. 
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FIG. 7: RerAD(t) (black line) and ImTAo{t) (gray line) with level broadening for the CS from the LUMO of naphthalene to 
the LUMO of TONE in PhCN. The level broadening is chosen as (a) T = 1 eV, (b) F = 0.1 eV, (c) T = 0.01 eV, (d) T = 0.001 
eV. 
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FIG. 8: RerAD(t) and ImTAD(t) with level broadening for the OS from the LUMO of naphthalene to the LUMO of TONE in 
PhCN. The gray line is the second-order perturbation result, and the black line is the full T-matrix result. The level broadening 
is chosen as (a) T = 1 eV, (b) T = 0.1 eV, (c) T = 0.01 eV, (d) T = 0.001 eV. 
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FIG. 9: Ci^{t) (black solid line), Ci^^it) (gray solid line), C^Uf) (black solid line), and Cf^"(t) (gray solid line) for the CS 
from the LUMO of naphthalene to the LUMO of TONE in PhCN. The black dashed lines are the approximation results for 
C'ae"(t) and Cf^'^it) using Eqs. ^ and The level broadening is chosen as (a) T = 1 eV, (b) T = 0.1 eV, (c) T = 0.01 

eV, (d) r = 0.001 eV. 
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FIG. 11: Time dependence of sin^ ^ (dark gray solid line) in comparison with (black dashed line) and |(0B|'!/'+)r 

(black dashed line), and cos^ 6 (gray solid line) in comparison with (black dotted line) and K^bIV'-)!^ (black dotted 

line). 
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FIG. 12: Energy gap dependence of the ET rate. The gray sohd, black sohd, and black dashed lines are the results for the 
quantum correction factors of Eqs. (|68[) . (|69p . and (|70[). respectively. The level broadening is chosen as (a) F = 1 eV, (b) 
F = 0.1 eV, (c) F = 0.01 eV, (d) F = 0.001 eV. 



